fl i rf  H u?> izuvc h  Reitilw:  o/ 1 E.tpitJhft* m  -  J-5,  JW/ 


MULTI-RESOLUTION  GLOBAL  MODELS  FOR  TELESEISMIC  AND  REGIONAL 

EVENT  LOCATION 

Michael  Antolik,  Goran  Ekstrom,  Adam  M.  Dziewonski,  Lapo  Boschi,  Bogdan  Kustowski,  and 

Jianfeng  Pan 

Harvard  University 

Sponsored  by  Defense  Threat  Reduction  Agency 
Contract  No.:  DSWAO 1-97- 1-00 17 


ABSTRACT 

Reconciliation  of  structures  shown  in  global  models  of  Earth  structure  with  those  shown  in 
detailed,  high-resolution  regional  models  has  become  one  of  the  most  important  tasks  in  seismic 
tomography.  Recent  studies  have  shown  that  combination  of  regional  seismic  phases  with 
teleseismic  phases  can  greatly  increase  location  accuracy;  however,  for  this  to  be  the  case,  the 
velocity  structure  has  to  be  well  resolved  on  a  fine  scale.  The  most  efficient  method  to  accomplish 
this  task  is  development  of  multi-resolution  models  using  local  (but  smooth)  basis  functions  that 
increase  in  resolution  from  regions  with  sparse  data  coverage  to  regions  of  dense  coverage.  This 
paper  reports  on  recent  results  from  an  ongoing  project  designed  to  bridge  the  gap  between  regional 
and  global  seismic  tomography  and  achieve  correspondingly  incremental  improvements  in  event 
location.  We  have  previously  developed  a  moderately  high-resolution  global  model  (equivalent  to 
spherical  harmonic  degree  18)  using  local  spline  functions.  This  model  improves  locations  using 
teleseismic  phases  over  other  more  finely  parameterized  models  (  Antolik  et  al.,  2001).  Here  we 
extend  the  same  techniques  to  simultaneous  development  of  more  detailed  regional  models  for 
inclusion  within  the  global  model.  We  describe  progress  in  constructing  a  new  regional  model  of 
the  Africa/Mediterranean  region  which  makes  use  of  surface  wave  dispersion  data,  regional  travel 
times  and  waveforms,  and  teleseismic  phase  arrivals.  This  model  also  incorporates  shear  wave 
anisotropy  and  agrees  well  with  other  published  models  for  the  region.  Other  regions  for  which  we 
intend  to  develop  combined  regional/global  models  include  the  former  Soviet  Union  and  North 
America. 

Tomographic  models  can  also  be  vastly  improved  through  the  use  of  sources  with  calibrated 
travel-time  data  (i.e.,  reference  events  with  accurately  known  locations).  However,  the  current 
database  of  such  events  still  suffers  from  sparse  coverage  (particularly  in  oceanic  regions).  A 
parallel  element  of  the  current  project  includes  improvement  of  event  locations  on  mid-ocean 
ridges  by  making  use  of  focal  mechanisms  and  accurate  bathymetry.  We  have  developed  a  database 
of  some  1500  large  events  and  are  expanding  to  smaller  events  using  the  Joint  Hypocentral 
Determination  technique.  We  present  the  results  along  with  experiments  designed  to  test  the 
accuracy  of  these  locations. 

KEY  WORDS:  Event  location;  seismic  tomography;  mantle  heterogeneity 
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OBJECTIVE 


The  aim  of  this  project  is  to  ultimately  improve  the  locations  of  earthquakes  and  other  seismically 
recorded  events  in  order  to  enhance  the  ability  to  monitor  the  Comprehensive 
Nuclear-Test-Ban-Treaty  (CTBT).  Our  strategy  is  based  on  developing  new,  detailed  3-D  models  of 
the  mantle,  with  an  emphasis  on  P  wave  structure.  This  involves  the  construction  of  global 
high-resolution  models  with  more  detail  in  certain  areas  where  particularly  good  data  coverage  is 
available.  In  our  approach,  3-D  velocity  structure  is  represented  in  terms  of  a  horizontally  layered 
and  laterally  varying  crust  of  variable  thickness,  overlying  a  mantle  whose  heterogeneous  P-  and 
S-velocity  structure  is  parameterized  radially  and  laterally  in  terms  of  cubic  B-spline  functions  of 
variable  density.  The  spline  functions  are  split  radially  at  the  670-km  discontinuity  in  order  to 
better  resolve  changes  in  heterogeneity  patterns  across  this  boundary  (Gu  et  al.,  2001).  The  shallow 
mantle  is  parameterized  as  an  anisotropic  medium  in  shear  velocity.  We  have  previously  developed 
a  new  joint  P  and  S  velocity  global  velocity  model  of  the  mantle  parameterized  in  spherical 
B-splines  (Antolik  et  al.,  2000).  This  model,  with  a  nominal  horizontal  resolution  equivalent  to 
degree  1 8  in  a  spherical  harmonic  expansion,  improves  the  average  mislocation  for  explosions  in 
Kennett  and  Engdhahl’s  (1991)  database  and  the  Prototype  International  Data  Centre’s  Reference 
Event  database  by  >20%  over  more  finely-parameterized  models  (see  Table  1).  We  are  now 
moving  to  take  advantage  of  better  data  coverage  obtainable  in  certain  regions  of  the  globe  by 
inverting  for  multi-resolution  models  with  finer  parameterization  confined  to  certain  regions.  The 
first  half  of  this  paper  describes  an  anisotropic  model  of  shear  wave  velocity  developed  for  the 
African/Mediterranean  region  developed  using  this  technique.  This  and  other  similar  models  will 
be  used  to  combine  both  regional  and  teleseismic  travel  times  for  the  improvement  of  event 
locations. 


Model 

Explosions 

Earthquakes 

PREM 

12.92 

18.83 

SP12 

7.83 

15.37 

MK12 

9.53 

17.40 

BDP98 

9.80 

17.13 

VWE97 

10.51 

17.33 

P362-17 

6.27 

13.75 

Table  1:  Average  mislocation  for  1 12  test  events  consisting  of  explosions  and  earthquakes  using 
PREM  (Dziewonski  and  Anderson,  1981)  and  five  3-D  mantle  models.  P362-17  is  our  new  global 
model  of  P  wave  velocity  described  in  Antolik  et  al.  (2000).  Only  teleseismic  P  wave  arrivals 
were  used  to  calculate  the  locations.  Location  errors  for  the  explosions  are  significantly  smaller 
than  for  the  earthquakes,  which  may  reflect  greater  errors  in  the  earthquake  reference  locations.  See 
Antolik  et  al.  (2001)  for  description  of  the  other  3-D  models  and  the  location  method  used. 

A  second,  subsidiary  objective  concerns  the  development  of  additional  techniques  used  in  locating 
events  teleseismically  and  regionally  with  sparse  datasets,  and  with  assessing  the  improvement  in 
accuracy  afforded  by  the  new  techniques  and  models.  To  address  this,  we  are  attempting  to 
improve  the  quality  of  locations  for  events  in  remote  regions  where  the  number  of  reference  events 
is  currently  small.  We  use  a  combination  constrained  inversion/JHD  technique  to  reconsider  the 
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locations  of  earthquakes  on  mid-ocean  ridges  and  transform  faults. 

RESEARCH  ACCOMPLISHED 


As  recent  tests  have  shown  that  compilation  of  accurate  travel  times  for  regional  (A  <  20°) 
first-arriving  phases  can  significantly  improve  location  accuracy  (Chen  and  Willemann,  2001),  we 
have  begun  to  make  rapid  progress  in  the  development  of  regional  velocity  models  of  both  P  and  S 
velocity,  including  consideration  of  anisotropy.  The  North  African/Mediterranean  region  is  one 
where  the  current  data  availability  is  such  that  velocities  in  the  upper  mantle  can  be  resolved  with 
higher  precision  than  that  obtainable  globally.  The  Mediterranean  Basin  is  a  geophysically 
complex  area,  governed  by  the  slow  collision  between  Africa  and  Eurasia.  The  boundary  between 
the  two  plates  is  not  entirely  well  defined.  Although  there  is  clear  geodetic  and  seismic  evidence  of 
northward  subduction  under  the  Aegean  and  Calabrian  arcs,  diffuse  seismicity  around  the  Adriatic 
Sea  and  Western  Mediterranean  requires  a  more  complex  description  of  plate  interaction,  which  is 
the  subject  of  current  research.  The  densely  populated  Italian  peninsula  has  a  long  history  of 
destructive  earthquakes,  and  a  better  understanding  of  the  regional  tectonics  would  have  a 
significant  role  in  the  context  of  seismic  hazard  mitigation  efforts.  Tomographic  studies  are 
naturally  a  means  to  this  end  (e.g.,  de  Jonge  et  al.,  1994;  Faccenna  et  al.,  2001). 

In  recent  years,  seismic  images  of  the  Mediterranean  Basin  have  been  derived  by  researchers  at 
University  of  Utrecht  (Wortel  and  Spakman,  2000)  and  at  I.N.G.V.  in  Italy  (Piromallo  and  Morelli, 
1997).  Both  groups  have  focused  on  explaining  observations  of  body- wave  (rather  than 
surface- wave)  travel  times;  the  Italian  group  have  emphasized  regional  data,  while  the  Dutch 
group,  with  the  recent  work  of  Bijwaard  (1999),  have  derived  a  high  resolution  regional  model 
within  the  framework  of  a  variable  resolution  model  encompassing  the  whole  globe. 


Figure  1:  Locations  of  horizontal  spline  knots  for  parameterization  of  the  multi-resolution  global  model  of 
anisotropic  shear  velocity. 
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Free  parameters  of  our  least-squares  inversions  are  lateral  variations  in  horizontally  and 
vertically-polarized  shear  velocities  throughout  the  upper  mantle,  described  as  linear  combinations 
of  radial  and  horizontal  cubic  splines  (Boschi,  2001).  The  basis  functions  are  more  densely 
distributed  within  the  Mediterranean  Basin,  to  achieve  a  higher  nominal  resolution.  Figure  1  shows 
the  horizontal  locations  of  the  spline  knots.  Globally,  the  resolution  of  the  inversion  approximates  a 
degree  18  spherical  harmonic  expansion  (as  in  the  inversion  for  joint  P  and  S  velocity  described 
above),  but  becomes  more  dense  within  the  region  of  interest.  The  measured  variations  in  phase 
delay  for  Rayleigh  and  Love  waves  between  periods  of  30  and  300  s  are  defined  as 

2  K  /*Aj  pa 

5$  j  (u)  =10^2^2  xik  /  Ki(r,Lo)fk(r,e,(j))drds,  (1) 

2  =  1  k  =  1  ,0  ,0 

where  the  index  j  represents  the  individual  measurements,  the  fk  represent  the  horizontal  and 
radial  spline  functions,  the  Xik  are  the  model  parameters,  and  the  index  i  refers  to  the  specific 
quantities  being  sought  in  the  inversion,  in  this  case  the  horizontal  and  vertical  shear  velocities.  The 
Ki  functions  are  known  as  the  “sensitivity  kernels”  (Anderson  and  Dziewonski,  1982)  or  partial 
derivatives.  An  original  feature  of  this  modeling  is  that  the  sensitivity  kernels  are  treated  as 
laterally  varying,  .i.e,  the  reference  model  used  is  not  one-dimensional.  Instead,  we  calculate  the 
partial  derivatives  for  model  CRUST5.1  (Mooney  et  al.,  1998)  overlying  a  PREM  upper  mantle. 
The  resulting  model  shows  significant  differences  in  absolute  amplitude  of  the  anomalies  when 
compared  to  models  constructed  in  the  conventional  manner  with  a  1-D  reference  model  (Boschi 
and  Ekstrom,  2001). 

Figure  2  shows  images  at  several  depths  of  the  global  model  of  vertical  shear  velocity  and  the 
regional  Mediterranean  model,  while  Figure  3  shows  a  cross-section  across  the  western 
Mediterranean  comparing  the  vsv  model  to  other  published  velocity  models.  Among  the  global 
features  shown  in  the  new  model  is  the  high-velocity  anomaly  in  vertical  shear  velocity  in  the  top 
of  of  the  upper  mantle  underneath  the  central  Pacific  first  reported  by  Ekstrom  and  Dziewonski 
(1998).  Our  new  Mediterranean  model  has  many  features  in  common  with  others  published  for  the 
region.  High-velocity  anomalies  associated  with  African  subduction  and  the  European  shield  can 
be  readily  seen  in  Figure  2,  as  well  as  low- velocity  models  associated  with  extension  (e.g.,  Aegean 
Sea).  The  cross  section  shown  in  Figure  3  indicates  that,  by  accounting  for  a  crust  of  laterally 
varying  properties  in  the  reference  model,  we  have  successfully  reproduced  the  strong  slow 
anomaly  (associated  with  tectonic  extension)  predicted  by  de  Jonge  et  al.  (1994);  in  the  same 
region,  where  the  crust  is  anomalously  thin,  Spakman  et  al.  (1993)  have  mapped  a  fictitious  fast 
anomaly.  The  thickness  of  the  crust,  whose  lateral  variations  Spakman  et  al.  (1993)  do  not 
consider,  trades  off  with  the  imaged  shallow  upper  mantle  velocities. 

We  are  currently  refining  our  imaging  techniques.  While  the  model  described  above  only  accounts 
for  measurements  of  surface-wave  phase  velocity,  we  expect  to  be  able  in  the  near  future  to  derive  a 
single  model  from  the  simultaneous  inversion  of  surface-wave  phase  velocity,  surface-wave  group 
velocity  (from  a  local  data  set  assembled  and  made  public  by  the  Lawrence  Livermore  National 
Laboratory),  body- wave  travel  times,  and  mantle  waveforms,  including  joint  consideration  of  P 
and  S  velocities.  This  step  should  improve  the  effective  resolution  of  our  image  in  the  upper 
mantle,  and  provide  a  reasonable  coverage  of  the  lower  mantle. 

In  addition  to  the  construction  of  more  detailed  and  accurate  global  velocity  models,  we  have 
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Figure  2:  Global  (left)  and  local  (right)  percent  heterogeneity  Svsv/vsv  at  six  depths  in  the  upper  mantle 
from  the  variable  resolution  model  of  anisotropic  shear  velocity. 
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Figure  3:  Comparison  between  four  different  images  of  the  upper  mantle,  shown  (in  cross-section) 
in  the  region  underlying  the  Western  Mediterranean;  the  left  panels  are  from  de  Jonge  et  al.  (1994). 
DRC-I  and  DRW  are  theoretical  P-velocity  models,  computed  by  de  Jonge  et  al.  (1994)  on  the 
basis  of  Dercourt  et  al.’s  (1986)  and  Dewey  et  al.’s  (1989)  tectonic  reconstructions,  respectively. 
EUR89B  is  the  tomographic  P-model  of  Spakman  et  al.  (1993),  based  on  regional  and  teleseismic 
P-travel  time  measurements,  and  including  a  simplified  1-D  reference  crustal  model.  On  the  right 
is  our  own  vsv  model  derived  from  surface  wave  phase  velocity  observations. 
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focused  our  efforts  on  building  a  larger  database  of  reference  or  “ground-truth”  events  to  be  used 
for  model  calibration.  So  far,  the  geographical  distribution  of  events  with  locations  known  to  an 
accuracy  of  5  km  or  better  is  extremely  limited.  Many  of  these  events  are  nuclear  explosions  and 
are  concentrated  in  only  a  few  source  regions,  making  comprehensive  testing  of  velocity  models 
and  calibration  of  new  seismic  stations  very  difficult.  We  have  established  a  catalog  of  ~  1500  large 
(M  >  5.5)  earthquakes  located  along  mid-ocean  ridge  and  transform  faults  using  a  constrained 
inversion  technique  where  the  event  is  confined  to  the  nearest  appropriate  plate  boundary  feature 
(ridge  or  transform  fault)  based  on  the  CMT  focal  mechanism.  The  inversion  problem  is  thus 
reduced  to  determination  of  a  single  distance.  Figure  4  shows  the  global  distribution  of  these  large 
events.  We  refer  to  these  events  as  “master  events”  in  the  discussion  that  follows. 


Figure  4:  Global  locations  of  relocated  master  events  from  the  CMT  catalog  along  mid-ocean 
ridges  and  transforms. 

Improvement  in  the  location  of  smaller  events,  for  which  a  focal  mechanism  has  not  been 
determined,  might  be  obtained  using  the  Joint  Hypocentral  Determination  (JHD)  technique.  To  do 
this,  we  form  clusters  of  events  using  those  events  in  the  ISC  catalog  within  a  circle  of  radius  100 
km  about  a  master  events.  The  differences  between  our  inversion  and  those  of  other  JHD 
procedures  (e.g.,  Pavlis  and  Booker,  1983;  Jordan  and  Sverdrup,  1981)  are  that  we  solve  for  source 
parameter  and  station  corrections  simultaneously,  and  that  the  locations  of  the  master  events  are 
held  fixed  during  the  inversion  process.  Stations  for  which  there  are  fewer  than  3  observations  for  a 
cluster  are  not  included,  and  we  use  the  LSQR  algorithm  (Paige  and  Saunders,  1982)  to  perform 
the  inversion  in  order  to  take  advantage  of  the  sparsity  of  the  inner  product  matrix.  The  depths  of 
all  the  events  are  assumed  shallow  and  are  fixed  at  10  km.  Figure  5  shows  the  results  of  one  JHD 
inversion  along  the  Romanche  Fracture  Zone  in  the  central  Atlantic.  After  inversion  the  events  are 
drawn  toward  the  fracture  zone  represented  by  the  low  in  the  bathymetry. 
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Additionally,  we  combined  the  inversion  for  each  master  event  together  to  solve  for  relocation 
vectors  and  station  corrections  for  a  regional  cluster  of  earthquakes  (i.e.,  JHD  with  multiple  master 
events).  Just  as  in  the  single-master-event  inversion,  each  master  event  influences  only  those  events 
within  a  radius  of  1 00  km.  In  the  combined  matrices,  we  apply  a  smoothing  weighting  function  to 
the  phases  of  each  master  event 

Wj  =J2e~10x(Ajk)2  ’  (2) 

K 

where  A ^  is  the  distance  between  the  jth  and  kth  master  events  in  the  region  under  study  and  the 
sum  is  over  all  of  the  K  master  events.  The  phase  weights  for  the  smaller  non-master  events  in  the 
inversion  are  fixed  at  unity,  and  all  events  for  which  fewer  than  30  P  phases  exist  are  omitted.  We 
again  use  an  iterative,  LSQR  algorithm  to  solve  the  equations.  This  process  has  been  applied  to 
approximately  1100  total  events  within  the  central  Atlantic  between  30°  S  and  30°  N. 


Figure  5:  Results  of  JHD  cluster  analysis  with  a  single  master  event  (open  diamond).  Lines  show 
relocation  vectors  of  the  smaller  events  in  the  cluster,  with  the  new  location  given  by  the  small 
diamonds.  The  vectors  have  been  enlarged  by  a  factor  of  two  in  length.  The  location  of  the  master 
event  was  held  fixed.  Bathymetry  is  taken  from  Smith  and  Sandwell  (1997). 


Figure  6  shows  the  regional  inversion  locations  for  148  events  along  the  Romanche  Fracture  Zone 
(40  master  events).  The  relocation  vectors  are  shown  in  the  bottom  panel,  and  indicate  that  most  of 
the  epicenters  become  better  aligned  along  the  linear  transform  fault,  except  near  the  eastern  edge 
where  there  is  considerable  complexity  in  the  intersection  of  the  transform  fault  and  adjoining 
ridge.  In  this  region  the  seismicity  shows  significant  deviation  from  the  great  circle  defining  the 
plate  boundary.  The  maximum  movement  obtained  for  the  smaller  events  in  the  JHD  process  (~50 
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km)  is  comparable  to  that  observed  for  the  original  master  events.  We  find  that,  after  the  JHD 
inversion,  the  majority  of  the  smaller  events  show  improvement  in  the  residual  RMS  (w.r.t.  the 
PREM  model)  over  the  initial  ISC  location  (the  ISC  uses  the  Jeffreys-Bullen  travel  time  tables). 
After  3  iterations,  a  reduction  in  the  weighted  RMS  residual  of  ~50%  is  achieved  for  the 
Romanche  Fracture  Zone  events.  The  significant  variance  reduction  achieved  for  the  residuals 
suggests  that  the  accuracy  of  the  new  locations  is  higher  than  the  original  catalog  locations.  The 
locations  obtained  in  this  step  can  be  compared  with  those  obtained  using  the  single-master-event 
cluster  analysis  in  order  to  assess  their  reliability. 
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Figure  6:  Regional,  combined  JHD  inversion  for  148  total  events  along  the  Romanche  Fracture 
Zone.  Top  panel  shows  the  initial  ISC  locations  for  the  smaller  events  in  the  inversion,  the  middle 
panel  the  new  locations,  and  the  bottom  panel  shows  the  relocation  vectors.  The  40  master  events 
included  (which  were  fixed)  are  not  shown. 
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CONCLUSIONS  AND  RECOMMENDATIONS 


Although  the  current  generation  of  global  velocity  models  provides  substantial  improvement  over 
1-D  velocity  models  for  teleseismic  earthquake  location  (Smith  and  Ekstrom,  1996;  Antolik  et  al., 
2001),  there  is  still  much  room  for  improvement.  The  consistent  accurate  location  of  small  events 
to  within  the  area  specified  for  on-site  inspection  under  the  CTBT  (1000  km2)  will  likely  require 
the  use  of  regional  phases  in  addition  to  teleseismic  travel  times,  and  also  possibly  consideration  of 
anisotropy  or  arrival  angle  data.  To  this  end,  we  are  developing  new  global  models  which  combine 
a  local  parameterization  with  the  advantages  of  smoothness  over  large  distances.  This  will  enable 
straightforward  replacement  of  portions  of  the  global  model  by  a  more  detailed  model  covering  a 
particular  region.  The  smooth  parameterization  will  allow  simple  computation  of  3-D  raypaths. 

Our  detailed  modeling  of  velocities  in  the  Mediterranean  shows  the  promise  of  simultaneous 
consideration  of  regional  and  teleseismic  data  in  the  tomographic  inversion  problem. 

Testing  of  new  global  and  regional  velocity  models,  as  well  as  calibration  of  newer  IMS  stations, 
will  benefit  from  the  compilation  of  reference  events  in  remote  regions.  Our  database  of  relocated 
events  on  mid-ocean  plate  boundaries  will  also  enable  construction  of  more  accurate  travel  time 
and  phase  velocity  datasets.  Since  published  locations  in  these  regions  are  known  to  contain  large 
errors,  this  database  may  be  of  considerable  use  even  if  the  relocated  epicenters  are  accurate  only  to 
within  1 0  km. 
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